# SLE human patients
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

# examine mva pathway ----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

# examine interested cytokines ------
key_cytokine <- c('Ccl20','Tslp','Flt3lg','Csf2','Tnf','Il6') |>
  str_to_upper()

vip.gene <- c(key_cytokine, 'HMGCS1', 'HMGCR', 'MVK')

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg') |>
  str_to_upper()

# import datasets -------------
## zheng2022 --------
zh22_dir <- list.dirs('mission/FPP/xiangya_sle_scRNA/data/zheng2022/',
                      full.names = T, recursive = F)

zh22.epi <- zh22_dir |> str_subset('E$')

zh22.epi

zh22.epi <- zh22.epi |>
  basename() |>
  set_names(zh22.epi, nm = _)

# 2 batches have different genes list, need use intersection
batch1_zh22 <- zh22.epi[1:4] |>
  Read10X(strip.suffix = T)

batch2_zh22 <- zh22.epi[5:6] |>
  Read10X(strip.suffix = T)

batch1_zh22 |> glimpse()
batch2_zh22 |> glimpse()

## copic hc epidermis -----
copic.path <- list.files('mission/FPP/xiangya_sle_scRNA/data/GSE275846_copic/',
                     full.names = T, recursive = F)

copic.mex <- copic.path |>
  read_geo_supp(name_regex = 'GSM\\d+')

copic.mex <- copic.mex |>
  GetAssayData()

## nakamizo hc epidermis --------
nkmz.path <- list.dirs('mission/FPP/xiangya_sle_scRNA/data/GSE281449_nakamizo/',
                         full.names = T, recursive = F)

### select matched samples -----
nkmz.path <- nkmz.path[c(2,3,5)]

nkmz.mex <- nkmz.path |>
  str_extract('age.+|Old.+|young.+') |>
  set_names(x = nkmz.path, nm = _) |>
  Read10X()

## merge --------
intsc_genes <- c(batch1_zh22, batch2_zh22, nkmz.mex) |>
  map(rownames) |>
  reduce(intersect)

intsc_genes |> glimpse()

lst.mex <- c(batch1_zh22, batch2_zh22, nkmz.mex) |>
  map(\(x)x[intsc_genes, ])

merged.mex <- lst.mex[[1]] |>
  RowMergeSparseMatrices(lst.mex[-1])

merged.mex |> glimpse()

sobj <- merged.mex |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(batch = case_when(str_detect(orig.ident, 'GSM') ~ 'copic',
                           str_detect(orig.ident, 'SLE\\dE') ~ 'zheng22_b2',
                           str_detect(orig.ident, 'SLE|HC') ~ 'zheng22_b1',
                           .default = 'nakamizu'),
         group = ifelse(str_detect(orig.ident, 'SLE'), 'SLE', 'HC'))

sobj |>
  dplyr::count(batch, group, orig.ident) |>
  DT::datatable()

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident', 'batch'))

sobj |>
  DimPlot(group.by = 'batch')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca.main')

sobj |>
  DimPlot(group.by = 'hpca.main', cols = 'Paired')

### determine V3h V3l KC --------------
sobj <- sobj |>
  mutate(manual.main = case_when(hpca.main == 'Neurons' ~ 'Melanocytes',
                                 hpca.main == 'Epithelial_cells' ~ 'Keratinocytes',
                                 .default = hpca.main))

sobj |>
  DimPlot(group.by = 'manual.main', cols = 'Paired')

sobj.kc4 <- sobj |>
  filter(str_detect(manual.main, 'Kera'))

sobj.kc4 <- sobj.kc4 |>
  quick_process_seurat(batch = c('orig.ident', 'batch'), skip_norm = T)

#### KC clusters umap ---------
sobj.kc4 |>
  DimPlot(cols = DiscretePalette(36),
          label = T, label.size = 2, repel = T, raster = T) +
  ggtitle('Leiden clusters of KC') +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/mix3.epi.sle.kc.cluster.umap.pdf')

sobj.kc4 |>
  filter(group == 'HC') |>
  mutate(seurat_clusters = as.character(seurat_clusters) |>
           fct_relevel('22','6','12','6')) |>
  DotPlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'), cols = 'RdYlBu') +
  theme_bw()

v3mk <- last_plot() |>
  pluck('data') 

v3mk |>
  mutate(id = fct_reorder(id, avg.exp)) |>
  ggplot(aes(features.plot, id, size = pct.exp)) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_pubr(legend = "right") +
  geom_point(aes(fill = avg.exp.scaled), shape = 21, stroke = .3) +
  labs(fill = "Average expression", size = "Percent expressed",
       x = 'Gene', y = 'Cell type')

sobj.kc4 |>
  filter(group == 'HC') |>
  FindAllMarkers(features = 'TRPV3', only.pos = T)

vip.fc <- sobj.kc4 |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'group',
                       features = vip.gene, ident.1 = 'SLE')

vip.fc |>
  filter(p_val_adj < .05) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape = 21) +
  scale_fill_gradient2(low = 'blue', high = 'red')

vip.fc |> DT::datatable()

sobj.kc4 <- sobj.kc4 |>
  mutate(v3.fine = ifelse(seurat_clusters %in% c(22,6,12), 'V3h', 'V3l'))

v3f.fc <- sobj.kc4 |>
  FindMarkersAcrossVar(split.by = 'v3.fine', group.by = 'group',
                       ident.1 = 'SLE',
                       features = c(kegg_mva, key_cytokine))

v3f.fc |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape = 21, stroke = .3) +
  labs(title = 'more v3h') +
  scale_fill_gradient2(low = 'blue', high = 'red')

sobj.kc4 |> write_rds('mission/FPP/xiangya_sle_scRNA/mix3.epi.kc.rds')

sobj.kc4 <- read_rds('mission/FPP/xiangya_sle_scRNA/mix3.epi.kc.rds')

sobj$v3.fine <- NULL

sobj <- sobj.kc4 |>
  as_tibble() |>
  select(.cell, v3.fine) |>
  left_join(x = sobj, y = _) |>
  mutate(manual.fine = ifelse(is.na(v3.fine), manual.main, v3.fine))

sobj <- sobj |>
  mutate(manual.fine = case_when(manual.fine == 'V3l' ~ 'TRPV3-lo-KC',
                                 manual.fine == 'V3h' ~ 'TRPV3-hi-KC',
                                 .default = manual.fine) |>
           fct_relevel(manual.fine, 'TRPV3-hi-KC', after = Inf),
         manual.main = fct_relevel(manual.main, 'Keratinocytes', after = Inf))

## main cell type umap -------------
sobj |>
  DimPlot(group.by = 'manual.fine',
          label = T, label.size = 2, repel = T, raster = T) +
  ggtitle('Cell types of epidermis') +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/mix3.epi.sle.celltype.umap.pdf')

# rds checkpoint ---------
sobj |> write_rds('mission/FPP/xiangya_sle_scRNA/mix3.epi.sle.rds')

sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/mix3.epi.sle.rds')

# figures -------------
## volcano: V3h KC SLE vs HC ----------
fine.type.slevhc.deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual.fine', group.by = 'group',
                       ident.1 = 'SLE')

fine.type.slevhc.deg |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/mix3.fine.slevhc.deg.csv')

fine.type.slevhc.deg <-
  read_csv('mission/FPP/xiangya_sle_scRNA/results/mix3.fine.slevhc.deg.csv')

fine.type.slevhc.deg |>
  filter(cluster == 'TRPV3-hi-KC', p_val_adj != 0 | gene %in% key_cytokine) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  plot_pub_volc(highlights = vip.gene,
                 group1 = 'SLE V3-hi-KC', group2 = 'HC V3-hi-KC', force = T)

publish_pdf('mission/FPP/figures/mix3.sle.epi.v3h.sle.vs.hc.pdf', width = 65)

## volcano: KC SLE vs HC ----------
kc.slevhc.deg <- sobj |>
  filter(str_detect(manual.fine, 'KC')) |>
  FindMarkers(group.by = 'group', ident.1 = 'SLE') |>
  as_tibble(rownames = 'gene')

kc.slevhc.deg |>
  filter(p_val_adj != 0 | gene %in% vip.gene) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  plot_pub_volc(highlights = vip.gene,
                group1 = 'SLE KC', group2 = 'HC KC', force = T)

publish_pdf('mission/FPP/figures/mix3.sle.epi.kc.sle.vs.hc.pdf', width = 65)

kc.slevhc.deg |>
  mutate(cluster = 'Keratinocytes') |>
  bind_rows(fine.type.slevhc.deg) |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/mix3.fine.slevhc.deg.csv')

## volcano SLE V3h vs V3l -----------
sle.v3hvl.deg <- sobj |>
  filter(group == 'SLE') |>
  FindMarkers(group.by = 'manual.fine', ident.1 = 'TRPV3-hi-KC',
              ident.2 = 'TRPV3-lo-KC') |>
  as_tibble(rownames = 'gene')

sle.v3hvl.deg |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/mix3.epi.sle.v3hvl.deg.csv')

sle.v3hvl.deg |>
  filter(p_val_adj != 0 | gene %in% vip.gene) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  plot_pub_volc(highlights = vip.gene, force = T,
                group1 = 'SLE TRPV3-hi-KC', group2 = 'SLE TRPV3-lo-KC')

publish_pdf('mission/FPP/figures/mix3.sle.epi.v3h.vs.v3l.volcano.pdf', width = 65)

## SLE MVA expr ------------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(kegg_mva, group.by = 'manual.fine') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.sle.epi.mva.expr.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.epi.mva.expr.csv')

### total KC ----------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(kegg_mva, group.by = 'manual.main') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.sle.epi.mva.total.kc.expr.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.epi.mva.expr.total.kc.csv')

## HC MVA expr ------------
sobj |>
  filter(group == 'HC') |>
  BubblePlot(kegg_mva, group.by = 'manual.fine') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.sle.epi.mva.expr.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.hc.epi.mva.expr.csv')

## SLE APF expr ----------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(key_cytokine, group.by = 'manual.fine') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis() +
  labs(title = 'SLE')

publish_pdf('mission/FPP/figures/mix3.sle.epi.apf.expr.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.epi.apf.expr.csv')

### total KC ----------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(key_cytokine, group.by = 'manual.main') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.sle.epi.apf.total.kc.expr.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.epi.apf.expr.total.kc.csv')

## HC APF expr ----------
sobj |>
  filter(group == 'HC') |>
  BubblePlot(key_cytokine, group.by = 'manual.fine') +
  theme_jpub +
  scale_size(range = c(0,4)) +
  RotatedAxis() +
  labs(title = 'HC')

publish_pdf('mission/FPP/figures/mix3.hc.epi.apf.expr.pdf', width = 60)

## KC differentiation marker expr ----------
sobj |>
  BubblePlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'), group.by = 'manual.fine') +
  theme_jpub +
  scale_size(range = c(0,3)) +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.sle.epi.KC.diff.marker.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  arrange(desc(id)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.epi.kc.diff.marker.expr.csv')

## MVA SLE vs HC logfc -------------
fine.type.slevhc.deg |>
  filter(!str_detect(cluster, 'V3')) |>
  inner_join(mva_fct) |>
  mutate(p_val_adj = if_else(p_val_adj == 0, 1e-300, p_val_adj),
         cluster = fct_relevel(cluster, 'Keratinocytes', after = Inf)) |>
  ggplot(aes(ordered, cluster, size = -log10(p_val_adj), fill = avg_log2FC)) +
  scale_fill_gradientn(colors = c('blue','white','red'),
                       values = c(0, 5, 1)) +
  theme_pubr(legend = "right") +
  labs(x = 'Gene', y = 'Cell type') +
  geom_point(shape = 21, stroke = .3) +
  scale_size(range = c(0,4)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.epi.slevshc.mva.logfc.total.kc.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  arrange(ordered) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.vs.hc.epi.mva.total.kc.csv')

## APF SLE vs HC logfc ---------
fine.type.slevhc.apf <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual.main', group.by = 'group',
                       ident.1 = 'SLE', features = key_cytokine,
                       logfc.threshold = 0, min.pct = 0)

fine.type.slevhc.apf |>
  mutate(p_val_adj = if_else(p_val_adj == 0, 1e-300, p_val_adj),
         cluster = fct_relevel(cluster, 'Keratinocytes', after = Inf)) |>
  ggplot(aes(gene, cluster, size = -log10(p_val_adj), fill = avg_log2FC)) +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(NULL,NULL)) +
  labs(x = 'Gene', y = 'Cell type') +
  geom_point(shape = 21, stroke = .3) +
  scale_size(range = c(0,4)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mix3.epi.slevshc.apf.logfc.total.kc.pdf', width = 60)

fine.type.slevhc.apf |>
  arrange(desc(cluster)) |>
  write_csv('mission/FPP/pub_source_data/mix3.sle.vs.hc.apf.apf.total.kc.csv')

kc.slevhc.apf <- sobj |>
  filter(str_detect(manual.fine, 'KC')) |>
  FindMarkers(group.by = 'group', ident.1 = 'SLE',
              features = key_cytokine, logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene')

kc.slevhc.apf |>
  mutate(p_val_adj = if_else(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes('SLE vs HC (Keratinocytes)', gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(NULL,NULL)) +
  labs(x = 'Gene', y = 'Cell type') +
  geom_point(shape = 21, stroke = .3) +
  scale_size(range = c(0,4)) +
  theme_jpub

publish_pdf('mission/FPP/figures/mix3.epi.kc.slevshc.apf.logfc.pdf')

## expr of V3 ----------
sobj.kc4 |>
  bill.violin('TRPV3', group.by = group)

sobj.kc4 |>
  get_abundance_sc_wide('TRPV3') |>
  left_join(x = sobj.kc4, y = _) |>
  ggplot(aes(group, TRPV3)) +
  geom_violin(aes(fill = group), scale = 'width') +
  stat_summary(fun = logtpm.mean, geom = 'crossbar',
               width = .3, color = 'black') +
  facet_wrap(~v3.fine, scales = 'free_y', ncol = 1) +
  labs(y = 'Normalized expression') +
  theme_jpub +
  scale_fill_manual(values = c('blue','red')) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3.7) +
  annotate(geom = 'text', x = 1.5, y = 3.7, label = '***') +
  expand_limits(y = 3.9)

publish_pdf('mission/FPP/figures/mix3.kc.v3.sle.hc.violin.pdf')

sobj.kc4 |>
  FindMarkersAcrossVar(split.by = 'v3.fine', group.by = 'group',
                       ident.1 = 'SLE', features = 'TRPV3')

## fraction of DC & KC --------
sample.conf <- sobj |>
  calc_frac_conf_on_grouped_count(orig.ident, manual.fine)

sample.conf |>
  filter(str_detect(manual.fine, 'KC|DC')) |>
  mutate(group = ifelse(str_detect(orig.ident, 'SLE'), 'SLE', 'HC'),
         sample = orig.ident, group.sum = NULL, conf.high = NULL,
         conf.low = NULL,
         .keep = 'unused')  |>
  write_csv('mission/FPP/pub_source_data/mix3.DC.KC.fraction.csv')

sample.conf |>
  filter(str_detect(manual.fine, 'KC|DC')) |>
  mutate(group = ifelse(str_detect(orig.ident, 'SLE'), 'SLE', 'HC')) |>
  ggplot(aes(group, fraction)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~manual.fine)
